Effect of fluctuations on the superfluid-supersolid phase transition on the lattice 
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We derive a controlled expansion into mean field plus fluctuations for the extended Bose-Hubbard 
model, involving interactions with many neighbors on an arbitrary periodic lattice, and study the 
superfluid-supersolid phase transition. Near the critical point, the impact of (thermal and quantum) 
fluctuations on top of the mean field grows, which entails striking effects, such as negative superfluid 
densities and thermodynamical instability of the superfluid phase - earlier as expected from mean- 
field dynamics. We also predict the existence of long-lived "supercooled" states with anomalously 
large quantum fluctuations. 
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I. INTRODUCTION 

The question of whether macroscopic quantum coher- 
ence can prevail in the presence of periodic order, ulti- 
mately leading to the existence of a supersolid, has been 
intensely debated since five decades [3, 0, S|. Of late, 
this topic has seen a renewed surge of interest, partly 
due to the observations in Q indicating potential sig- 
natures of a supersolid phase of 4 He. Because of the 
inherent complexity of 4 He, it is useful to gain further 
understanding by studying supersolid phases in other sys- 
tems - such as the Bose-Hubbard model, which can be 
realized experimentally via cold bosonic atoms in opti- 
cal lattices @. For on-site interactions only, the phase 
diagram at T = contains the superfluid and the Mott 
insulator state @, 0|. Adding interactions across sites 
(next nearest or higher) in a so-called extended Bose- 
Hubbard model, a superfluid-supersolid phase transition 
may occur [1, @, [13, EH, [12] in addition to further Mott- 
like phases. Here the term supersolid is associated to 
an order parameter (with a well-defined phase) which is, 
in contrast to the homogeneous superfluid ground state, 
not the same for all lattice sites, but periodically modu- 
lated. At the heart of supersolid formation is the generic 
phenomenon that an instability towards density modu- 
lations occurs if the excitation spectrum dips below zero 
for a finite wavevector. 

Within mean-field theory, i.e., neglecting all fluctua- 
tions, properties of the supersolid phase were studied 
in [13j. However, the evanescent excitation energies at 
the transition suggest that (thermal and quantum) fluc- 
tuations should play an important role near the critical 
point. The impact of these fluctuations can be taken into 
account with quantum Monte Carlo simulations, see, e.g., 
@,[l3|- Despite the strength of this method, these simula- 
tions are always restricted to a specific (low-dimensional) 
lattice of finite size and a small sample of the full Hilbert 
space. In the following, we consider an arbitrary peri- 
odic lattice and develop an analytic expansion into mean 
field plus fluctuations where the size of the fluctuations 



and the validity of the expansion is controlled by a small 
parameter. Therefore, our derivation is complementary 
to other numerical and analytical approaches using for 
example duality to vortex field theory [13]. To the end 
of devising a controlled mean field expansion, we begin 
by introducing the concept of weighted operator sums in 
the next section. 



II. WEIGHTED OPERATOR SUMS 

We consider the extended Bose-Hubbard model on an 
arbitrary lattice as described by the Hamiltonian 

H = ^ [ T <*P "1^/3 + 7^"/3 a^a^aaaA , (1) 

a/3 \ 1 J 

where a, [3 label the lattice sites and a) a ,a@ are the as- 
sociated bosonic creation/annihilation operators. The 
kinetic term is determined by the hopping matrix T a p 
and the interaction part by V a p . (Both matrices are real 
and symmetric.) Since the above Hamiltonian cannot be 
diagonalized analytically, we have to employ some ap- 
proximations. To this end, we introduce the concept of 
weighted operator sums defined via 

*s[f] = t|t E Uala a ), (2) 

with a set S C N of \S\ elements a £ S and a func- 
tion f a of order one, i.e., which does not scale with \S\. 
Hence the limit limigi^oo Xs[f] exists (in an appropri- 
ate sense, e.g., as a weak limit) while all single addends 
are suppressed by 1/|5| for large \S\. Examples for the 
form |(2]) include all (local) one-site operators such as 
Xr Q \[l] = a a for |5| = 1 as well as the (global) Fourier 
components ^2 a a a exp{ika}/L = ak/VTj with |5| = L 
being the total number of sites for a one-dimensional 
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chain S = [1,L]. Now, considering the commutator be- 
tween two such weighted operator sums 



X s [f],Xs>[f] 



\sns'\ 
\s\ x \S'\ 



X 



sns- 



(3) 



with fa(at,,a a ) = [f a (a a , a a ), a a )], we find that 

they are suppressed for large \S\ due to \S (1 S'\ < 
min{\S\, \S'\}. Hence the limit ]im|g|_ >0Q .X' 1 s[/] com- 
mutes with all other weighted operator sums (including 
all local operators) and can thus be approximated by a 
c-number within the relevant Hilbert space generated by 
weighted operator sums acting on the ground (or ther- 
mal) state. This motivates the following asymptotic ex- 
pansion for large \S\ ^> 1 



X s [f} = C [f} + 



]S] \ s \ 



(4) 



where the leading term Co[/] can be approximated by 
a c-number and the sub-leading operators C'i/2[/] and 
Ci/ 2 [f] generate the commutator J3]), of order 1/|5|. 

Applying the concept of weighted operator sums to 
operators like = J2^ap/L ° r other Fourier com- 

ponents, we arrive at the mean-field expansion 



(5) 



where ip a denotes the mean field and corresponds to 
the leading parts Cq[J] in Eq. (j4]) while the fluctua- 
tions x a with (\a) — incorporate the non-commuting 
remainders. Note that (in contrast to [lH) the filling 
n a = (a^cia) = |0^| + (XaXa) is here not assumed to be 
large; \ip a \ ls the condensate part and (x a Xa) is the re- 
maining thermal or quantum depletion. Hence the fluc- 
tuations Xa are not necessarily small compared to the 
mean field tp a ; e.g., for half-filling n a = 1/2, the vari- 
ance is obviously of order one. In order to simplify the 
full equation of motion derived from (TJ [h = 1] 



id t a a = (T a/3 ap + V a phpa a ) 



(6) 



P 



we assume that the interaction V a p involves a large num- 
ber D ^> 1 of sites P on a roughly equal footing. This 
could be the case, for example, for long-range interac- 
tions or for a large number of spatial dimensions. For 
normalized potentials ~Y^gV a p = Vs = 0(1), we may 
then apply the concept of weighted operator sums (J2j) to 
the term ^ a V a php and obtain 



V <*P h P 

P 



P 



0(1/VD) 



(7) 



from Eq. (J4j) - However, one must be careful: simply re- 
placing hp by np in ((6]), we would lose the phonon modes. 
The sub-leading term 0(1/VD) can only be neglected if 
there is no other small (or large) term involved. This is 



precisely the case for modes with long wavelengths over 
many lattice sites, where the sum over T a pap, for ex- 
ample, is also very small and hence the 0(1/VD) con- 
tributions become relevant. In order to describe long- 
wavelength modes correctly, we insert Eq. J5|) into Eq. 
to obtain the Gross-Pitaevski! equation 

id t ^ a = (r<xpi>p + v a3 [iVvsl 2 + (xlxp}] 0«) , (8) 



where we have replaced J2 g VapXpXp by its expectation 
value according to the above arguments, plus the remain- 
ing fluctuation part 



dtXa = J2(T a pX0 + V a p \ipp\ 2 + (xpXp 



+ V C 



a/3 [i>*pX0 + i>0Xp (0a + Xa) 



(9) 



Again, the second line is suppressed by 0(1/ vD) and 
will only be relevant for long-wavelength modes, which 
involve a sum over many sites a. In this case, however, 
the c-number term J2 a ipa will dominate the fluctuation 
term J2 a Xa in view of Eq. (@1) and hence we may ap- 
proximate the bracket (0 a + Xa) in the second line by 
ip a , arriving at a linear operator equation 



idtXa = X! [TapXP + V a/3 |0/?| 2 + {x\x,P 



Xa 



+ V«P [rpXP + <P0XI] 0a) + 0(1/ VD) , (10) 

which corresponds to the Bogoliubov-de Gennes equa- 
tions for the fluctuations. Note that the approximation 
from Eq. |(9|) to Eq. lfT7J|) neglects the exchange of par- 
ticles between the condensate \ip a \ and the thermal or 
quantum depletion (x a Xa) ■ This exchange is governed 
by the sub-leading term J2^V a ptp^(xpXa) , which could 
be added to Eq. (JHJ - 



III. QUASIPARTICLE MODES 

In order to introduce quasiparticle modes, we assume 
translational invariance, i.e., that T a p and V a p only de- 
pend on the distance a — (3 and that the condensate den- 
sity is homogeneous |0 a | — 101- Nevertheless, we may 
still have a constant phase gradient t] in our sample, i.e., 
we set ip a — |0| exp{— ifit + irja}. In this case, we may 
diagonalize Eq. (fl0|) via a Fourier transformation 



idtXk = (T k+ri + V s n + V k \i/j 2 \) X k + V k ip 2 X l k ■ (11) 

Note that in more than one spatial dimension, a and (3 
as well as k and r\ will be multi-indices (labeling the real 
and the inverse lattice, respectively) and t}ol is a scalar 
product. Assuming reflection invariance T k = T k = T_ k 
and Vfc = V k * = V~ k for the lattice, we see that this 
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symmetry fc — > — A; is broken for the modes Xk by the 
phase gradient 77. The quasiparticle Hamiltonian 

Hx = Y,{xK Tk +v + V v n +\^\ 2 Vk)xk 

k ^ 

+Y (V' 2 xix t -, + h.c.)} = E^ I (i2) 



the condensate fraction \ip 2 \ and thus weakens the term 
If/^lVfe in Eq. lfl4|) responsible for the roton dip. Ergo, 
heating up the supersolid state may yield the superfluid 
phase (as long as the condensate does not disappear al- 
together ip = 0), which will become important for the 
discussion of "supercooled" states we turn to in section 

El 



can be diagonalized via the Bogoliubov transformation 
Xk = u k b k + VkbL k with \uf.\ — \v%\ = 1. This yields the 
Bogoliubov coefficients 



1 



h 



iik 



Vk 



l-*r 



Ik = \ wi + 2w k -1-Wk, w k 



VkW 



(13) 



For Wk — —2, the coefficients diverge due to Ik = 1 (lead- 
ing to the instability for 77 = 0, to be discussed below). 
The quasiparticle frequencies obey the dispersion relation 

(T k=0 = 0) 



-(T k+V - r fc _) ± JT 2 + 2\^V k T k , (14) 



where T k = (Tk+r) + Tk- v )/2 and thus the branches are 
connected by u k = —w~k- 

In the continuum limit, i.e., for small k <C 1, we may 
approximate Tk ~ k 2 /(2m) due to Tk = T k — TL^ and 
Tfc = o = with the mass m being determined by the hop- 
ping rates. For small phase gradients 77 <C k -C 1, we 
then reproduce the usual Galilei shift 



(ut+vk) 2 = \Tp 2 \V k — + 



(2m) 



(15) 



where v = 77/771 is superfluid velocity. Now, even for 
purely positive V a /3, the Fourier transform Vk may be- 
come negative for some k and hence the dispersion re- 
lation may develop dips (similar to the roton dip in su- 
perfluid 4 He) . If Vk is sufficiently negative (compared to 
Tk), the dispersion curve u>k may even dive below zero. 
Ignoring the fluctuations discussed below, the onset of in- 
stability, LOk = 0, marks the end of the (homogeneous) su- 
perfluid phase and the beginning of the supersolid phase 
where \ip a \ is periodic, i.e., inhomogeneous. The phase 
gradient 77 favors the supersolid phase, i.e., the transition 
superfluid — * supersolid occurs earlier for non-vanishing 
77. For 77 = T = 0, the frequencies ujk=± k , at the roton 
wavenumber become imaginary beyond the critical point 
and hence these modes start to grow exponentially. For 
77 > and T = 0, the transition occurs earlier and and is 
slower since the frequency u> k =+k, becomes negative, but 
not imaginary. Hence only the coupling to some environ- 
ment (fixed by the lattice) induces an instability of these 
quasiparticle modes. On the other hand, the depletion 
(XaXa) due to thermal or quantum fluctuations favors 
the superfluid phase since it reduces (for a fixed filling n) 



IV. SUPERFLUID DENSITY 

Now let us study the response of the system to a small 
phase gradient a a — > a Q exp{i77a}, which determines the 
superfluid fraction. The interaction part \V a p a a a^a a a,j3 
of the Hamiltonian ([I]) does not change, but the kinetic 
term yields 



' a/3 



(16) 



which is a measure for the total current dH/drj oc J, 
cf. the Fourier expansion in (|12p. E.g., for T a p oc 
S a ,(3+i + S a ,j3-i — 25 a ,i3 [HI, we get the usual expres- 



sion J oc i(a\ 



h.c). In the continuum limit of 



T fe sa k /(2m), we obtain 



J oc 



dH 
dr] 



1 

m 



E 



t\W" 



■XkXk 



J 4 , + J x , (17) 



where the first term rj\ip 2 \ in the bracket is the condensate 
(i.e., mean-field) contribution and the second one, J x , 
stems from the fluctuations. Inserting the Bogoliubov 
transformation Xk — Ukbk + Vk$_ k , we find 



f k + \^ 2 \Vk 
2Jf 2 + 2f k \iP 2 \V k 



(xl fc X-fc)o, (18) 



i.e., the expectation value of the fluctuation part in the 
ground state vanishes (J x )o — 0. Even though the quasi- 
particle frequencies are different in opposite directions 
for a non- vanishing phase gradient 77, Uk 7^ oj-k, the Bo- 
goliubov coefficients are still symmetric \uk\ = \u-k\ and 
\vk\ = \v-k\. Because of the symmetry Ik = l~k, from 
Eq. lfl3|), quantum depletion does not contribute to the 
current [171 ] . 

In a thermal ensemble, as described by the density ma- 
trix g = e~xjp{—H x /T}/Z, however, quasiparticle modes 
with Uk 7^ u)-k will have different occupation numbers 
and hence we do get a contribution to the total flux from 
the fluctuations (J) oc X)fc( r ?l' ( /' 2 I Clearly, 
near the superfluid-supersolid phase transition, the ma- 
jor contributions occur around the roton minima at ±fc*. 
Here, we consider for simplicity one spatial dimension 
only, but the main results apply to higher dimensions as 
well. Let us first study the case 77 = 0. In the continuum 
limit k <C 1, we may use a Taylor expansion 



uj 2 = 2T k \^ 2 \Vk + T 2 



+ j 2 (k - K 



(19) 
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around the roton minimum at the critical wavenumber 
fe* 1, where 7 is the curvature of the roton dip. Ap- 
proaching the phase transition corresponds to the limit 
lo 2 — > and for small with <C T, the leading term 



scales as 



L ^ 

k 



O 



Tina;,,, 
7 



(20) 



At the critical point lu* — 0, the fc-integral over the ther- 
mal distribution (i>\b k ) ps T/w^ becomes weakly diver- 
gent near the roton dip at where u>k ps 7|fc — fe*|, 
leading to the logarithmic singularity lnw*. 

Now, adding a small phase gradient, one roton mini- 
mum is lifted and the other one approaches the u = 
axis even closer. Hence the thermal quasiparticle occu- 
pation numbers (b\.bk) react in opposite ways and induce 
a net current - which is opposite to the condensate flux 
^l^ 2 !- The change due to — > u>* ± vk* scales as 



(Jx 



TyK 

w*7 



(21) 



For small enough w» or, alternatively, for large enough 
temperatures T > T cr = <D(uj*mj\?p 2 \/kl), the cur- 
rent induced by the thermal fluctuations can easily com- 
pensate the condensate (mean-field) contribution r)\ij} \, 
Thus, the superfluid fraction can be significantly reduced 
- and may even become negative (which is also occurring 
in 7r-Josephson junctions [la]), i.e., the phase gradient 77 
entails a net current (J) in the opposite direction. 

Such a negative superfluid density induces a thermo- 
dynamical instability As discussed before, the cur- 
rent (J) is a measure for the response of the system to 
a phase gradient (dH/drj). Inserting the canonical en- 
semble q = exp{—H/T}/Z, we see that the expectation 
value (dH/drj) = dF/drj equals the change of the free 
energy F = E — TS = (H) + T(lng). Since a stable equi- 
librium state in an isothermal environment corresponds 
to a minimum of the free energy, a negative superfluid 
density d(J)/drj < shows that the system is unstable 
against the spontaneous generation of local phase gradi- 
ents (since rj = is a maximum of the free energy) . 

Note that a negative superfluid density does not re- 
quire a large thermal depletion: as we may infer from 
Eq. (|20| . the thermal occupation number ^2 k (pJ>k) scales 
merely logarithmically with uj* and hence it can be much 
smaller than the condensate fraction \ip 2 \ (e.g., for T -C 7 
and fc* <C 1 fc*lnaj* < 1). Of course, in addi- 
tion to thermal occupation, the condensate is also de- 
pleted by quantum effects. This quantum depletion sur- 
vives at zero temperatures and is given by Eq. ifTSj) via 
(XaXa) = J2k(xlxk)/L. With the same approximations 
as in Eq. ([20|) . we get again merely a logarithmic depen- 
dence 



(xlx a )o = O 



k 2 In a;* 

rwy 



(22) 



Consequently, a vanishing superfluid density d(J)/dr) = 
0(t]), which marks the end of the (homogeneous) super- 
fluid phase may occur even when the total (thermal plus 
quantum) depletion is very small \ip\ 2 ^> (XpXp) an d 
hence the condensate fraction is still near one. Note that 
this behavior is opposite to (bulk) superfluid 4 He, where 
the superfluid fraction (near 100% for low temperatures) 
strongly exceeds the condensate part (of order 10%). 



V. "SUPERCOOLED" STATES 

Although the depletion was small \ip\ 2 ^> (xpXp) m the 
cases under consideration, we would like to stress that the 
presented controlled mean-field expansion can also be 

applied to the case of large depletions (xpXp) — ^(IV 7 ! 2 )- 
This generalization can be achieved by demanding that 
Vk is strongly peaked at the origin Vk=o — Ve — 0(1) 
and much smaller otherwise Vjfei>k D <C 1 such that the 
width kg of the peak at the origin is much smaller than 
the typical k- values (position and breadth 1/^/7) as- 
sociated with the roton-dips (where (x k Xk) yields the 
major contribution). To see how this works, let us 
compare J2p VafiXpj which must be small within our 
approach, with the depletion (xlx.a) = Y,k(XkXk)/L. 
Calculating the squared norm (| Va/3Xp\ 2 ), we get 

J2k \Vk\(XkXk)/L. Similarly, higher orders yield a sum 
over several wavenumbers containing Fourier components 
at linear combinations of roton wave-numbers Vk±k' ■ 
Consequently, all these terms are suppressed even though 
the quantum depletion may be large. 

Given these requirements, one may obtain "super- 
cooled" states, which are long-lived superfluid phases in 
a parameter region where the true ground state is su- 
persolid. In order to demonstrate the main idea, let us 
consider the following gedanken experiment: We start in 
the superfluid phase at T — 0, where 90% of the parti- 
cles are in the condensate |i/' 2 | and 10% in the quantum 
depletion (XaXa)- Now we remove 80% of the particles 
(e.g., by a Raman transition with no momentum transfer 
incurred) by decreasing the condensate part \ip 2 \ only, 
i.e., we leave the modes with k 7^ forming the quan- 
tum depletion untouched. Simultaneously, we increase 
the interaction strength Vk (e.g., via a Feshbach reso- 
nance) such that the product |V ,2 |^4 remains constant, 
leaving the quasiparticle spectrum intact. After that pro- 
cedure, half of the remaining particles are in the conden- 
sate \ip 2 \ and the other half are in the quantum depletion 
(XaXa) = IV' 2 !- These anomalously large quantum fluctu- 
ations are caused by the increased interaction Vk, which 
is so strong that the true ground state (with this filling 
n = \ip 2 \ + (XaXa)), having significantly smaller deple- 
tion, is supersolid. However, the immediate transition to 
the supersolid state is prevented by the fact that only 
half the particles are in the condensate. Because the 
quasiparticle modes have the same positive energies as 
before, the system is linearly stable. Similar to the ther- 
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modynamical instability caused by a negative superfluid 
density, the decay to the true supersolid ground state is 
mediated by the sub-dominant term Yip VapippixpXa) ■ 
Ergo, the predicted "supercooled" state is long-lived and 
thus might be accessible to an experimental verification. 

VI. CONCLUSION 

In summary, by means of a controlled expansion into 
powers of the small parameter 1 / \f~D, yielding the mean 
field ip phis (thermal and quantum) fluctuations Xa, we 
are able to study the impact of these fluctuations onto 
the superfluid-supersolid phase transition analytically. In 
addition to the instabilities indicating the end of the 
(homogeneous) superfluid phase known from mean field 
dynamics, which occur when the roton dip touches the 
lo = axis, the fluctuations induce a thermodynamic in- 
stability even before reaching the classical critical point 
cj» =0. This breakdown of the homogeneous superfluid is 



associated with a negative superfluid density and occurs 
rather slowly, since changes of the mean field %p induced 
by fluctuations Xa are governed by the sub-dominant 
term Yip Va/3^(X/3Xa) in Eq. (fTOf . which is effectively 
a 0(l/-\/Z))-correction to the Gross-Pitaevskii Eq. ([8]). 
Finally, even though the thermodynamical instability ef- 
fect is governed by thermal fluctuations, quantum fluc- 
tuations do also generate intriguing phenomena near the 
critical point like supercooled states. 
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